####################################################################################
# Ajustes dos modelos multinomiais a partir dos dados do Enade
# Autoria: Leonardo Rodrigues
# Código escrito e executado no RStudio (Version 1.3.1073)
# Tempo aproximado para rodar o script: 30 minutos
####################################################################################

start_time <- Sys.time()

### 1 Preparação ----
# Leitura dos pacotes
library(foreign)
library(nnet)
library(stargazer)
library(EMT)
library(DescTools)
library(broom)
library(pscl)
library(mlogit)
library(reshape2)
library(dplyr)
library(stringr)
library(ggplot2)
library(forcats)
library(ggpubr)
library(gt)
library(Amelia)
library(tidyr)

#Leitura do banco
enade <- read.table("data/originaisEnade/enade.txt", sep="\t")

#Ajustes das variáveis
table(enade$cor) #Categoria "Outros" será excluída N: 41885.

#Mudando nome das categorias para aproveitar script anterior
enade$cor = ifelse(enade$cor == "Brancos", "Branca",
                   ifelse(enade$cor == "Negros", "Não-Branca", NA))

#renomeando variáveis para aproveitar script anterior
#estabelecendo a categoria de referência para a modelagem
enade$year <- as.factor(enade$ciclo)
enade$year = relevel(enade$year, ref = "2")
enade$tese = ifelse(str_detect(enade$curso,"ENGENHARIA"), enade$curso, "outros")
enade11 <- enade
enade11$year <- enade11$ciclo
enade11$year <- as.factor(enade11$year)
enade11$year = relevel(enade11$year, ref = "2")
enade11$grupo2 <- enade11$diploma
enade11$grupo2 <- as.factor(enade11$grupo2)
enade11$grupo2 = relevel(enade11$grupo2, ref = "Engenharias")
enade11$regiao <- enade11$regiao1
enade <- NULL 


#### 2 Obtendo os ajustes dos modelos ao acrescentar cada variável ----


# Primeiro modelo: área ~ idade
rm(enade)
max_iterations <- 500

#criando data.frame para salvar os ajustes
ajuste <- data.frame("Tipo de Ajuste" = c("Log-Lik Intercept", "Deviance", "R2 McFadden",
                                          "AIC", "BIC", "Log-Lik", "R2 McFadden Adj", "DF"),
                     "Idade" = NA, "I+R" = NA, "I+R+S" = NA,
                     "I+R+S+R" = NA, "I+R+S+R+E" = NA, "I+R+S+R+E*A" = NA)

m1 <- multinom(grupo2 ~ age, data = enade11, maxit = max_iterations)

ajuste$Idade[1] <- PseudoR2(m1, "logLik0")
ajuste$Idade[2] <- m1$deviance
ajuste$Idade[3] <- PseudoR2(m1, "McFadden")
ajuste$Idade[4] <- m1$AIC
ajuste$Idade[5] <- BIC(m1)
ajuste$Idade[6] <- nnet:::logLik.multinom(m1)
ajuste$Idade[7] <- PseudoR2(m1, "McFaddenAdj")
ajuste$Idade[8] <- m1$edf

rm(m1)

#Segundo modelo: área ~ idade + região
m2 <- multinom(grupo2 ~ age + regiao, data = enade11, maxit = max_iterations)

ajuste$I.R[1] <- PseudoR2(m2, "logLik0")
ajuste$I.R[2] <- m2$deviance
ajuste$I.R[3] <- PseudoR2(m2, "McFadden")
ajuste$I.R[4] <- m2$AIC
ajuste$I.R[5] <- BIC(m2)
ajuste$I.R[6] <- nnet:::logLik.multinom(m2)
ajuste$I.R[7] <- PseudoR2(m2, "McFaddenAdj")
ajuste$I.R[8] <- m2$edf

rm(m2)

#Terceiro modelo: área ~ idade + região + sexo
m3 <- multinom(grupo2 ~ age + regiao + sexo, data = enade11, maxit = max_iterations)

ajuste$I.R.S[1] <- PseudoR2(m3, "logLik0")
ajuste$I.R.S[2] <- m3$deviance
ajuste$I.R.S[3] <- PseudoR2(m3, "McFadden")
ajuste$I.R.S[4] <- m3$AIC
ajuste$I.R.S[5] <- BIC(m3)
ajuste$I.R.S[6] <- nnet:::logLik.multinom(m3)
ajuste$I.R.S[7] <- PseudoR2(m3, "McFaddenAdj")
ajuste$I.R.S[8] <- m3$edf
rm(m3)

#Quarto modelo: área ~ idade + região + sexo + cor
m4 <- multinom(grupo2 ~ age + regiao + sexo + cor, data = enade11, maxit = max_iterations)

ajuste$I.R.S.R[1] <- PseudoR2(m4, "logLik0")
ajuste$I.R.S.R[2] <- m4$deviance
ajuste$I.R.S.R[3] <- PseudoR2(m4, "McFadden")
ajuste$I.R.S.R[4] <- m4$AIC
ajuste$I.R.S.R[5] <- BIC(m4)
ajuste$I.R.S.R[6] <- nnet:::logLik.multinom(m4)
ajuste$I.R.S.R[7] <- PseudoR2(m4, "McFaddenAdj")
ajuste$I.R.S.R[8] <- m4$edf
rm(m4)

#Quinto modelo: área ~ idade + região + sexo + cor + escolaridade dos pais
m5 <- multinom(grupo2 ~ age + regiao + sexo + cor + edu, data = enade11, maxit = max_iterations)

ajuste$I.R.S.R.E[1] <- PseudoR2(m5, "logLik0")
ajuste$I.R.S.R.E[2] <- m5$deviance
ajuste$I.R.S.R.E[3] <- PseudoR2(m5, "McFadden")
ajuste$I.R.S.R.E[4] <- m5$AIC
ajuste$I.R.S.R.E[5] <- BIC(m5)
ajuste$I.R.S.R.E[6] <- nnet:::logLik.multinom(m5)
ajuste$I.R.S.R.E[7] <- PseudoR2(m5, "McFaddenAdj")
ajuste$I.R.S.R.E[8] <- m5$edf
rm(m5)

#Sexto modelo: área ~ idade + região + sexo + cor + escolaridade dos pais e interação com ciclo
m6 <- multinom(grupo2 ~ age + regiao + sexo*year + cor*year + edu*year, data = enade11, maxit = max_iterations)

ajuste$I.R.S.R.E.A[1] <- PseudoR2(m6, "logLik0")
ajuste$I.R.S.R.E.A[2] <- m6$deviance
ajuste$I.R.S.R.E.A[3] <- PseudoR2(m6, "McFadden")
ajuste$I.R.S.R.E.A[4] <- m6$AIC
ajuste$I.R.S.R.E.A[5] <- BIC(m6)
ajuste$I.R.S.R.E.A[6] <- nnet:::logLik.multinom(m6)
ajuste$I.R.S.R.E.A[7] <- PseudoR2(m6, "McFaddenAdj")
ajuste$I.R.S.R.E.A[8] <- m6$edf
rm(m6)

#salvando os resultados
write.table(ajuste, "data/output/ajusteareas.txt", sep="\t")


### 3 Criando tabela dos ajustes do modelo ---

ajustes <- read.table("data/output/ajusteareas.txt", sep="\t")

tabela <- ajustes %>%
  gt(rowname_col = "Tipo.de.Ajuste") %>%
  tab_header(title = md("Estatísticas de Ajuste do Modelo")) %>%
  cols_align(align = "center") %>%
  fmt_number(vars("Idade","I.R","I.R.S","I.R.S.R","I.R.S.R.E","I.R.S.R.E.A"),
             rows = c(1),
             decimals = 3, use_seps = FALSE, dec_mark = ",") %>%
  fmt_number(vars("Idade","I.R","I.R.S","I.R.S.R","I.R.S.R.E","I.R.S.R.E.A"),
             rows = c(2,3,4,5,6,7),
             decimals = 5, use_seps = FALSE, dec_mark = ",")

end_time <- Sys.time()

end_time - start_time